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Abstract. In this work it is shown how the immersed boundary method of (|64T) for modeling 
flexible structures immersed in a fluid can be extended to include thermal fluctuations. A stochastic 
numerical method is proposed which deals with stiffness in the system of equations by handling 
systematically the statistical contributions of the fastest dynamics of the fluid and immersed struc- 
tures over long time steps. An important feature of the numerical method is that time steps can 
be taken in which the degrees of freedom of the fluid are completely underresolved, partially re- 
solved, or fully resolved while retaining a good level of accuracy. Error estimates in each of these 
regimes are given for the method. A number of theoretical and numerical checks are furthermore 
performed to assess its physical fidelity. For a conservative force, the method is found to simulate 
particles with the correct Boltzmann equilibrium statistics. It is shown in three dimensions that 
the diffusion of immersed particles simulated with the method has the correct scaling in the physical 
parameters. The method is also shown to reproduce a well-known hydrodynamic effect of a Brownian 
particle in which the velocity autocorrelation function exhibits an algebraic (t~ 3 / 2 ) decay for long 
times f(|; Il6t l2Ct |2a ; |37t l38t |54| : |67l ; |78|) . A few preliminary results are presented for more complex 
systems which demonstrate some potential application areas of the method. 
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1. Introduction. In modeling many biological systems it is important to take 
into account the interaction of flexible structures with a fluid. The immersed boundary 
method of |64|) has found wide use as an efficient numerical method for simulating such 



systems. Some examples include the study of blood flow around heart valves (65), 

Hi- 



wave propagation in the inner ear (|32f). and the generation of lift in insect flight 
With experimental advances in molecular and cellular biology has come an increasing 
interest in developing methods to model qualitatively and quantitatively microscopic 
biological processes at the cellular and subcellular level (fl3 ; r27l l4lh . The immersed 
boundary method provides a promising framework for simulating such systems. 

At the cellular level the fluid may consist of either the aqueous environment 
outside of the cell or the cytoplasm within. Some important flexible structures in the 
cellular context include the outer cell membrane, intracellular vesicles, cytoskeletal 
fibers, and molecular motor proteins. These structures play an important role in cell 
motility or cell division among other processes H). 

In such systems the relevant features can span a range of length scales from 
tens of microns or more for the outer cell membrane and cytoskeletal fibers to tens 
of nanometers for individual cytoskeletal monomers and motor proteins. At these 
length scales thermal fluctuations of the system become significant and in many cases 
appear crucial to achieve biological function. Some examples include for ce g eneration 
and progression of molecular motors alo ng cyto skeletal fibers ([]]; |4l|; |84 ). osmotic 
effects such as vesicle and gel swelling (2(| 33; 51; 61; 83), and polymerization effects 
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involved in force generating processes in cell motility ( 34 ; 6(1 72; 81). 

Modeling such complex cellular systems in molecular detail is infeasible with 
methods such as molecular dynamics as a consequence of the immense computational 
cost required to resolve the broad range of active length and time scales. This suggests 
that a coarse-grained numerical approach must be taken which does not resolve all of 
the detailed physics but rather attempts to make approximations that yield effective 
equations to capture the most relevant features of the dynamical phenomenon being 
studied. Here we discuss how the framework of the immersed boundary method can 
be extended for use in such modeling by including thermal fluctuations to capture 
dynamical phenomena at the cellular length scale. 

The theory of nonequilibrium statistical mechanics indicates that the influence of 
thermal fluctuations on a mechanical system can typically be represented through the 
addition of thermal forcing terms which decorrelate rapidly in time. The forcing can 
then be represented by appropriate "white noise" processes. This generally involves 
a nontrivial structure of correlations between the state variables in such a way that 
there is an energy balance between the thermal forcing and dissipation of the system 
so that a corresponding fluctuation-dissipation theorem for the system is satisfied 
|I3;[ZQ|). 

Several computational fluid dynamical schemes have been extended toward the 
microscale through such an inclusion of thermal forces. The most widely used ap- 
proach is known as Stokesian or Brownian dynamics (fl5c 122c FBI l77h . In this approach 
the structures are modeled as collections of rigid "elementary particles" which inter- 
act through force laws derived by approximating the fluid dynamics by a quasi-steady 
Stokes flow. This latter approximation is strictly appropriate only when, among other 
assumptions, the fluid density is much less than the density of the structures (I19h . 
a condition better met in engineering applications (such as suspensions ( 14c 1 15c 1771)) 
than in physiological settings (fl3T) ■ The result of the underlying approximations in the 
Stokesian/Brownian dynamics method is a rather strongly coupled system of stochas- 
tic equations for the motion of the elementary particles. For a well-designed com- 
putation the cost of a simulation can be rendered roughly proportional to NlogN, 
where N is the number of elementary particles f77l). In the presence of fast time 
scales arising from thermal fluctuations and possibly chemically activated processes, 
the impact of the quasi-steady Stokes approximation and the representation of the 
elementary particles as rigid (rather than flexible) on the accuracy of the simulation 
is not yet clear (|73l ). 

Another approach for modelin g flu ids with immersed structures is Dissipative 
Particle Dynamics <S El H H ^SBl H3 IH HH ■ The method is built phenomeno- 



logically in terms of "fluid particles" which represent a parcel of fluid along with its 
collection of immersed structures. When thermal forces are included in the method 
the fluid particles arc simulated with a stochastic system of equations modeling their 
(soft) interactions. This method however does not readily extend to the microscopic 
domain since the immersed structures within a parcel are not resolved in detail. Dis- 
sipative Particle Dynamics may however be appropriate for somewhat larger scale 
simulations in which one is interested in the effects of a numerous collection of im- 
mersed polymers or other structures on the dynamics of a fluid flow. 

A different class of approaches which emphasize the role of the fluid dynam- 
ics while making other simplifications has also been proposed. These include finite- 
element ([761 ) and lattice-Boltzmann (|48l ) methods, in which the computational fluid 
dynamics are extended to include thermal forces in the fluid equations following the 
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framework of The immersed boundary method which we shall discuss belongs to 
this broad class of methods (|45T). A theoretical approach with certain similarities to 
the immersed boundary method with thermal fluctuations has been proposed in (62), 
but differs in how thermal forces are treated for the immersed structures. A virtue of 
the immersed boundary method when compared to other methods is the straightfor- 
ward physical manner in which it approximates the interaction of the fluid with the 
flexible structures. 

A key feature of the immersed boundary method, distinguishing it from Stoke- 
sian/Brownian dynamics, is that the dynamics of the fluid are represented in the 
immersed boundary equations so that subtle inertial effects of the fluid can be in- 
corporated into the thermally fluctuating dynamics. For example, as demonstrated 
in Subsection 15.21 the method captures the slow decay (t~ 3 / 2 ) in the tail of the au- 
tocorrelation function of the velocity of an immersed particle. Another advantage 
of tracking the fluid dynamics is the natural way in which the immersed boundary 
method can respect the topology of flexible structures so that, for example, polymers 
do not cross themselves or each other. This feature gives the immersed boundary 
method the potential for efficient simulation of polymer links and knots. 

A basic description of how thermal fluctuations can be incorporated within the 
immersed boundary method was presented in (|45h. and theoretical analysis of the 
physical behavior of the method through an asymptotic stochastic mode reduction 
calculation was developed in |46l). The immersed boundary method was found in 
these theoretical works to produce generally the correct physical behavior for the 
thermal fluctuations of immersed structures. 

Here we present a derivation for the thermal fluctuations of the immersed bound- 
ary method in the context of the time dependent Stokes equations. We then present 
a new numerical method developed from a novel time discretization of the stochastic 
equations. In addition, further theoretical analysis of the framework is performed to 
investigate the physical behavior of the method and comparisons are made between 
theory and numerical simulations. 

For stochastic differential equations most traditional finite difference methods, 
such as Runge-Kutta, achieve a lower order of accuracy than for deterministic ordinary 
differential equations as a consequence of the nondiffcrcntiability of Brownian motion 
and its order i 1 / 2 scaling in time fiil ). When considering the full system of equations 
of the immersed boundary method, these issues are further compounded by a wide 
range of length and time scales that arise in many problems. 

For small length scale systems in which the Reynolds number is small and the 
fluid flow is Stokesian to a good approximation, the time scales associated with the 
fine-scale fluid modes can be considerably faster than the time scales of the large- 
scale fluid modes and the immersed structures. In many problems it is the dynamics 
of the immersed structures and the large-scale features of the fluid flow that are of 
interest. The fine-scale features of the fluid are incorporated in simulations primarily 
to determine their effects on the larger scales but are often not of direct interest in 
and of themselves. 

Since the full system is rather stiff as a result of the fast time scales of the fine- 
scale fluid modes, we would like to be able to take time steps which underresolve those 
scales of the fluid which are not of primary interest in the simulation. However, we 
must take care in how this is done, because otherwise the effects of the underresolved 
fluid modes on the more interesting degrees of freedom can be misrepresented. For 
example, the simple use of a method such as Euler-Marayama (|44l ) leads to poor 
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accuracy for the stochastic dynamics when long time steps are used. In this case, the 
particle trajectories are found to be overly diffuse in the sense that the wrong scaling 
is obtained for the mean squared distance traveled. 

After a presentation of the general framework of the immersed boundary method 
without spatiotemporal discretization in Section^ a stochastic analysis is developed 
in Section [3] to design a numerical method which maintains accuracy over time steps 
which can be longer than the time scales of some or all of the fluid modes. The sys- 
tematically derived numerical scheme presented here improves upon in several ways 



a correction-factor approach taken in (|45l ) to achieve long time steps. The numeri 



cal method is constructed from a new time discretization of the immersed boundary 
equations, in which the equations are integrated analytically using standard tech- 
niques from stochastic calculus under well controlled approximations. The numerical 
method allows for the statistical contributions of the fast stochastic dynamics of the 
fluid, which are not explicitly resolved over long time steps, to be accounted for in 
the dynamics of the immersed structures. Moreover, the correlations in the statistics 
between the fluctuations of the degrees of freedom of the system are handled sys- 
tematically, allowing for consistent realizations of the velocity field of the fluid and 
immersed structures to be simulated. 

In Section [4j error estimates are given which indicate that the method attains a 
good level of accuracy (in a strong statistical sense (0)) whether the fluid modes are 
completely underresolved, partially resolved, or fully resolved. Only the degrees of 
freedom of the immersed structures constrain the time step. The numerical method 
handles a broad range of time steps in a unified manner, so that depending on the 
application, the fast dynamics of the fluid can either be explicitly resolved or under- 
resolved, with their effects correctly represented on the structural degrees of freedom. 

In Section [5l an expression for the diffusion coefficient of immersed particles is 
derived. The predictions are compared with the results of numerical simulations 
showing good agreement with the theory for different particle sizes and both short 
and long time steps. It is further shown that for intermediate time steps, the method 
captures a well-known hydrodynamic effect of a Brownian particle, in which the decay 
of the autocorrelation function of the velocity of the particle decays algebraically 
(t- 3/2 ) (§ [H; [H [H; [H; [H [H 111)- Also in Section El numerical simulations 



are performed which confirm that the method produces the correct osmotic pressure 
and equilibrium statistical distribution for the position of particles within an external 
potential at finite temperature. To demonstrate more complex applications of the 
method, simulations are then presented which investigate behavior of the osmotic 
pressure associated with confinement of molecular dimers and polymer knots in a 
microscopic chamber, as well as a basic model of a molecular motor protein immersed 
in a fluid subject to a hydrodynamic load force. 

The physical consistencies we demonstrate in the method through our analysis 
and numerical experiments indicate that the stochastic immersed boundary method 
is a viable means to model on a coarse scale the influence of thermal fluctuations on 
the interaction of fluids and flexible structures. This suggests that the method holds 
promise as an effective approach in modeling complex biological phenomena which 
operate at the cellular and subcellular level. 

2. Fluid- Particle Equations. For the physical systems with which we shall be 
concerned, the relevant length scales will typically be on the order of tens of microns 
or smaller. The amplitude of the velocity fluctuations on these scales are sufficiently 
small, relative to the viscosity and length scale, that the Reynolds number is very 
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small. This will allow us to neglect the nonlinear advection term in the Navier-Stokes 
equation for the fluid dynamics. However, we will not drop the time derivative term 
(as is often done in low Reynolds number limits (fHoh ) because the dynamical time 
scales arising from Brownian motion and possibly certain vibrational modes of the 
immersed structures are in general too fast to allow this. That is, the dynamics can 
generally exhibit time scales which are much shorter than the advection time scale 
(length scale divided by velocity scale), so we can drop the nonlinear advection term 
but not the time derivative of the velocity. This leads us to the time-dependent Stokes 
equations for an incompressible fluid, which read 

(2-1) p 9u g' l) = /iAu(x, t) - Vp + ftotalfo *) 

(2.2) Vu = 0, 

where p is the pressure arising from the incompressibility constraint, p is the fluid 
density, p, is the dynamic viscosity, and ftotal IS ^ ne total force density acting on the 
fluid. 

The force density acting on the fluid arises from two sources. The first source is the 
forces applied to the fluid by the immersed structures and particles. This component 
of the force density, denoted f pr t, generally arises from the elastic deformations of 
immersed structures, but they may also be applied externally and transmitted by the 
immersed structures to the fluid. The second contribution to the force density is from 
the thermal fluctuations of the system and is denoted by fthm- Each of these force 
densities will be discussed in greater detail below. Together, they comprise the total 
force density acting on the fluid 

(2-3) f total( x > f ) = f prt( x > t) + f t hm( x > *)• 

The immersed boundary model for fluid-structure and fluid-particle coupling 
treats the flexible structures and particles to first approximation as part of the fluid, 
representing their structural properties (such as elasticity) through the force density 



term f pr t (|64f). All structures, such as membranes, polymers, and particles, are 
modeled as a collection of M discrete "elementary particles," with locations denoted 
by {X^'l (t)}jL 1: which interact with force laws appropriate to their structural prop- 
erties. For simplicity, we consider the case in which all forces can be described in 
terms of a conservative potential V({X}) depending on the positions of the collection 
of elemen tary particles. More general force relations, including active forces, can be 



included (|64f). For notational convenience in the exposition, the range of indices for 
the collection of elementary particles will often be omitted. We will often refer to the 
forces exerted by the immersed structures as "particle forces," since the structures 
are represented in the numerical method as a collections of interacting particles. 

In the immersed boundary method, elementary particles of size a are represented 
by a function S a (x) which may be thought of as a Dirac delta function smoothed 
over a length scale a in such a manner that the smoothed delta function has good 
numerical properties (See Appendix |A1 and (f&j)). 

This smoothed delta function is used both in converting the force associated with 
an elementary particle to a localized force density acting on the fluid: 

M 

(2.4) fprt(x,t) = £(- V xB']V)({X(t)}) i a (x-X^(t)) 
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and in computing the velocity of the elementary particle by an interpolation of the 
fluid velocity in its vicinity: 

dXtr] (+) r 

(2.5) ~~dt = J 5 a (x - (t))u(x, t)dx. 

The integration is over the entire domain A of the fluid. In the immersed boundary 
formulation, this includes the space occupied by the immersed particles and structures, 
which are thought of as parts of the fluid in which additional forces happen to be 
applied. In particular, the domain A is independent of time, despite the motion of 
the immersed material. 

We also note that in the present context the parameter a is a physical parameter 
of the model, since it is supposed to represent a physical dimension of an elementary 
particle. In particular, a is not a numerical parameter which is supposed to vanish 
along with the meshwidth for the fluid computations as it is refined. In this respect, 
the use of smoothed delta functions described here is different from the standard use 
of such functions in immersed boundary computations. The idea that smoothed delta 
functions could be used to model the physical dimensions of immersed objects was 
previously proposed in (fHHh under the name "force cloud method". 

3. Numerical Method. We shall now discuss a numerical discretization and 
specification of the thermal force density for the equations 12.1112.51 defining the im- 
mersed boundary method. A summary of the numerical method in algorithmic form 
is given in Subsection 13. li followed by a heuristic discussion in Subsection 13.21 and a 
mathematical derivation in Subsection 13.31 

3.1. Summary of the Numerical Method. The numerical method is based 
upon a finite difference discretization of the differential equations 12.11 12.51 describing 
the coupled dynamics of the fluid and the immersed structures. The fluid variables 
(velocity field u and pressure field p) are represented on a periodic grid with length L 
along each direction, N grid points along each direction, and grid spacing Ax = L/N . 
The values of these fields on the lattice will be denoted through subscripted variables 
such as u m and p m , where the subscript m = (mi, m2, 7713) is a vector with integer 
components indicating the grid point in question (relative to some arbitrarily specified 
origin). The position of the grid point with index m is denoted by x m . 

The Discrete Fourier Transform (DFT) of the fluid variables plays an important 
role in the numerical simulation scheme, and is related to the physical space values 
on the grid through the formulas: 

(3.1) u k 

(3.2) u m 

where each of the sums in the above equations runs over the N 3 lattice points 
defined by < k^) < N — 1, and < < N — 1, where the parenthesized 

superscripts I = 1,2,3 denote the Cartesian components of the indicated vector. In 
fact any translate of these blocks of lattice points could be used equivalently in the 
sums, due to the underlying periodicity. 

Time is discretized into time steps At, and the values of the system variables at 
the n th time step, corresponding to the time t n — nAt, are denoted with a superscript 



= -^3 ^2 Um exp M 2?rk ' m/iV) 

m 

= Uk exp (i27rk • m/N) , 

k 
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integer n. The procedure by which these variables are updated from one time step to 
the next is now described: 

1. The structural forces exerted by the immersed structures are computed and 
the lattice values of the particle force density field which is applied to the 
fluid is obtained from 

N 

(3-3) fm = E - ( V x(^) ({X n }) 5 a (x m - X"^). 

j'=i 

Here and afterwards, we drop the subscript "prt" from the particle force 
density. The Fourier coefficients f k of this particle force density field are 
computed using a Fast Fourier Transform (FFT). 

2. The Fourier coefficients of the velocity field of the fluid are updated by the 
stochastic recurrence 



(3.4) u^ 1 = e-** At K + (1 - e- QkAt ) p^ + p£B k , 

where p k denotes the projection orthogonal to gk defined by 

(3.5) gjf } = sin(27rk« /N) /Ax 



which is used to enforce the incompressibility constraint 12.21 The factor 
H k = cr^k accounts for the thermal fluctuations over the time step, where r) k 
denotes a complex vector-valued random variable independent in k, having 
independent real and imaginary components, each of which are Gaussian 
random variables with mean zero and variance one. The variance of H k is 
determined in Subsection 13.3.41 and is given by 

(3.6) al = ^ (i _ cxp (_2o k At)) , 
where 

(3.7) «k = ^ - cos(27rk«/jV))) 

pAx 

and 



(3.8) Ac = 



with 

(3.9) /C={k | k w = or k (j) = N/2, j = 1, 2, 3| . 
3. The elementary particle positions are updated by 

(3.10) x n+1 'M - X n '[ J l = J2 ^( x m - X n '^)r™ Aa; 3 , 
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where is the time integrated velocity field of the fluid. It is obtained by 
a discrete Inverse Fast Fourier Transform (IFFT) of appropriately generated 
random variables T£ in Fourier space: 

(3.11) f u m (s)ds = Vf|J-exp(i27rk-m/iV). 

Jtn k 

The rjj are computed from 

(3.12) fjj = H k + ci, k ^C + c2,kPkG k) 

* 71 

where H k is obtained from step 2 and Hk is computed from steps 1 and 2 by 
(3.13) 

H k = i^l^MuC + (-+(-) 2 (ex P (-a k At) - 1)) P ~^. 

yak \akj J 

The random variable G k is computed from scratch for each mode k by gener- 
ating a complex vector- valued random variable having independent real and 
imaginary components, each of which are Gaussian random variables with 
mean zero and variance one. The constants in 13. 121 are given by 



(3.14) ci k = — tanh 

ak 

and 



1 ' 




(3.15) c 2 ,k = \l[ ^ ] ( a k At-2tanh^^^ 

In this manner the time integrated velocity field is consistently generated 
with the correct correlations with {u£} and {u£ +1 } from steps 1 and 2. For 
more details, see Subsections 13.3.51 and 13.3.61 
The computational complexity of the method, when excluding the application specific 
forces acting on the immersed structures, is dominated by the FFT and IFFT, which 
for a three dimensional lattice requires 

0(iV 3 log(iV)) arithmetic steps. 

3.2. Heuristic Discussion of the Numerical Method. We now briefly dis- 
cuss each step of the numerical scheme to give some intuition into how the method 
operates. A more rigorous mathematical discussion and derivation is given in Subsec- 
tion [231 

The first step of the numerical scheme computes the structural forces exerted 
by the elementary particles as a function of their configuration, and computes the 
discrete Fourier transform of the corresponding force density field acting on the fluid. 
The second step updates the fluid velocity in Fourier space by integrating over the 
structural forces and thermal forces experienced over a time step. The appearance 
of the time step At in exponential factors is due to the design of the method to 
maintain accuracy even if the fluid dynamics are partially resolved or underresolved 
by the time step, as in "exponential time differencing" schemes (HSlai S3). The key 



time scale paired against the time step in these formulas is l/ct k , which describes the 
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time scale of viscous damping of the Fourier mode k of the fluid, as simulated by the 
numerical method. The first term of the stochastic recurrence equation [33] represents 
the viscously dissipated contribution of the fluid velocity from the previous time step. 
The second term represents the contribution from the structural forces during the time 
step. Since the elementary particle dynamics are assumed to be resolved by the time 
step, the structural force density itself appears in a simple way as effectively constant 
over the time step. The response of the fluid velocity to this force has an exponential 
dependence on the time step to account for the possible levels of resolution of the 
viscous damping. The third term accounts for the thermal fluctuations over the time 
step through a mean zero Gaussian random variable H k with variance cr£ describing 
the magnitude of the net contribution of the thermal fluctuations to the fluid velocity 
over the time step, see equation 13.61 Note that for time steps At long compared to 
the relaxation time l/«k of the velocity Fourier mode, the variance approaches the 
constant value ksT/(pL 3 ) corresponding to the equilibrium equipartition value. On 
the other hand, for At <C 1/ak, the variance of the thermal velocity increment is 
proportional to At, so the magnitude of the increment is proportional to %/ At. This 
latter scaling is typical for the response of physical systems to noise driven by a large 
number of weak inputs (i.e., molecular fluctuations) (J44j). In both the second and third 
terms, the projection enforces the incompressibility of the fluid. The distinction 
in the definition of the factor with respect to wavenumbers k, as specified by 
the set /C, is a purely technical issue related to the discrete Fourier transform; see 
Subsection 13.3.31 

The third step of the numerical method updates the positions of the elementary 
particle positions composing the immersed structures. The random variable rep- 
resents the fluid velocity at lattice point x m integrated over the time step, and is 
generated in Fourier space using a procedure which also ensures that the immersed 
structures move with the correct correlations with the previously computed fluid ve- 
locity values uJJ and u£ +1 - The formulas defining arise from an exact formula 
for integrating the fluid velocity over a time step, under the assumption that the 
structural forces can be treated as constant over the time step. 

3.3. Derivation of the Numerical Method. The derivation first considers a 
spatial discretization of equation 12.11 while leaving the system of equations continuous 
in time to avoid technical issues associated with the continuum formulation of the 
stochastic immersed boundary method with thermal forcing (|46l ). Since the equations 
are meant to serve as a physical model for the dynamics of the immersed structures and 
fluid, a thermal forcing is derived for the semi-discretized system which is consistent 
with equilibrium statistical mechanics in Subsection 13.3.31 The time discretization of 
the numerical method for both the dynamics of the fluid and immersed structures is 
then discussed in Subsection 13.3.41 and Subsection 13.3.51 The method takes special 
care to account for correlations between the dynamics of the fluid and immersed 
structures, which is discussed in Subsection 13.3.61 We remark that throughout the 
derivation, the integration steps of the numerical method are designed to maintain 
accuracy even when the time step does not fully resolve the dynamics of the fluid 
modes. 

3.3.1. Semi-discretization. The equations 12.11 and 12.21 can be discretized in 
space by finite difference approximations for the spatial derivatives ((641 ) 



(3.16) 




9=1 
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Pm+ee Pm-ei M) , s 

2Ax +t total^ Xm '^ 



3 U W (t)-ll W (t) 

1 ' ' 2Ax 

where denotes the standard basis vector with all zero entries except for a one in 
the I th position. The parenthesized superscripts denote the vector component. 
The equations for the fluid-particle coupling become 

M 

(3.18) fprt(x mj i) = -(V x u'iV0({X(t)}) <5 a (x m - Xb"l(t)) 

i'=i 

(3.19) ^=U(X^),t) 

(3.20) U(x,t)=£ <5 a (x m - x)u(x m ,i)A 



x 3 . 



Since we do not take the limit a — > as the meshwidth is refined (sec above), 
U(X^l(t),t) does not become the same as u(XW(t),t), even in the limit Ax — > 0. 
That is, we are not simply evaluating the fluid velocity at X^(t), but instead av- 
eraging it over a region of width determined by the parameter a. This averaging 
procedure ensures that the particle velocity remains finite and well-defined in the con- 
tinuum limit, in which the pointwise values of the thermally fluctuating fluid velocity 
diverge. Indeed, according to general statistical mechanical principles for continuum 
fields |47]; 0; For ) (and not from any feature particular to the numerical method) , the 
thermally fluctuating fluid velocity field manifests increasingly wilder fluctuations on 
smaller scales, and in the theoretical continuum limit approaches a sort of white noise 
structure. This is not inherently problematic for physical interpretation, which only 
requires that meaningful values be obtained from averages over finite volumes, such 
as the size of a probe or an immersed structure. Such averages are indeed finite, as 
can be understood intuitively through central limit theorem considerations by viewing 
them as averages of a large number of independent mean zero random variables due 
to the rapid spatial decorrelation of the noisy continuum velocity field. From a more 
mathematical standpoint, an application of the convolution theorem to the definition 
of U in (|3.18D yields that U can be represented as a nicely convergent Fourier series 
due to the decay of the Fourier coefficients of the smooth function S a . 

To obtain other desirable behaviors, we also remark that it is important that we 
use the same weight function 5 a in averaging the fluid velocity as we do in applying 
force to the fluid, since this ensures that energy is properly conserved in the fluid- 
particle interaction. With appropriate care in the construction of 5 a , one can further 



ensure that momentum and angular momentum are conserved as well; see (|64l ) 



3.3.2. Fluid Equations in Fourier Space. The Stokes equation is given in 
Fourier space by 

(3.21) ^ = -a k u k - ip^Pkgk + p^ftotal.k 

(3.22) gk-u k = 0, 
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where 

(3.23) a k = Y (1 - cos(27rk(^ /JV))) 

P Ax ]~i 

(3.24) g[ j) = sin(27rk w /A0/Ax. 

Since the velocity field of the fluid is real-valued, a further condition that must 
be satisfied by solutions of the equations 13.211 - 13.221 is 



(3.25) u N -k = Uk, 

where N is shorthand for (N, N, N) T and the overbar denotes complex conjugation. 
Provided the force is real-valued, it can be shown that if this constraint holds for the 
initial conditions it will be satisfied for all time. 

The Fourier coefficients pk(t) of the pressure need to be chosen in order to ensure 
that the incompressibility constraint is satisfied. They can be determined by pro- 
jecting both sides of equation 13.211 onto gk. By the incompressibility constraint 13. 22[ 
both of the terms involving Uk and its time derivative are zero under the projection. 
This gives at each time 

,„ oa \ - u\ ^ k ' f total,k( < ) 

(3.26) Pk (t) = — — . 

For those values of k that make gk = 0, the incompressibility constraint is trivial, 
and by convention we shall take p k {t) = for such k. 

For future reference, let the projection in the direction gk be denoted by 

(3.27) A - |f 
and the projection orthogonal to gk be denoted by 
(3-28) p ^ = fj-^EL 

For those modes for which gk = 0, the corresponding projections will be understood 
to be defined pj[ = and pj^ = The set of indices on which gk = is given by 

(3.29) £ = {k | k< j) = or = N/2, j = 1, 2, 3} . 

For a function w m defined over the discrete lattice sites indexed by m, the cor- 
responding projection operations in physical space are given by 

(3.30) (p"w) m = ^ pj[wk exp (i27rk • m/JV) 

k 

and 

(3.31) (p^w)™ = ^2 Pk *k exp (i27rk • m/N) . 
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3.3.3. Thermal Forcing. Following standard practice in nonequilibrium sta- 
tistical mechanics (47; l70l ). the thermal fluctuations of the system are modeled as 



Gaussian white noise. Formally, the Fourier coefficients of the thermal forcing can be 
written as 



(3.32) f thmjk = p^2D^ 



dB k (t) 

dt 



The factor Z? k (which is to be specified) describes the strength of the thermal forcing of 
the k*' 1 mode and B k (i) denotes a complex- valued Brownian motion with the real and 
imaginaryparts of each component consisting of an independent standard Brownian 
motion (|44T ). The dependence on k will be discussed below. 

Standard Brownian motion B(t) for our purposes will refer to the continuous 
stochastic process which is defined by the following properties: 

(i) B(0) = 0, 

(ii) E{B(t 2 )-B(t 1 ))=Q, 

(iii) ^ (|S(i 2 ) - ^(ii)! 2 ) - |i 2 - 

(iv) The increments Bfo) — B(t\) and B(t^) — B(t%) are 
independent Gaussian random variables whenever t\ < t2 < £3 < £4, 

where E (•) denotes the expected value. Standard Brownian motion in d dimensions 
is defined as a stochastic process where each vector component is an independent one- 
dimensional Brownian motion. For a further discussion of the properties of Brownian 
motion and related technical issues, see (|59l ) or (|3ll ). 

The discretized Stokes equation 13.211 with only thermal forcing (no immersed 
structural forces) can be expressed in stochastic differential notation as 



(3.33) du k = [~a k u k - ip-^kgk] dt + v/2~Dk"dB k (t) 

(3.34) gk • u k = 

(3.35) UN-k = u k , 

where c?B k (t) denotes increments of the complex- valued Brownian motion associated 
with the k*' 1 mode. To ensure that the thermal forcing be real-valued, the Brownian 
increments are correlated in k by the constraint 

(3.36) c?B N _ k = dB k . 

As discussed in Section 13.3.21 the pressure can be expressed in terms of the force 
acting on the fluid using [3~26l By formal substitution into 13.331 the incompressibility 
constraint can be incorporated through an appropriate projection operation which 
allows for the two equations 13.331 and 13.341 to be expressed as the single equation 

(3.37) du k + a w u k dt = ^/^D^p^dB^t). 

Since the incompressibility constraint is equivalent to p k u k = u k , the constraint 
will be satisfied for all time provided it holds at the initial time. Consequently, when 
g k 7^ (k ^ JC as defined in !3.29|) . the real and imaginary part of the stochastic process 
u k (£) remain in the plane orthogonal to g k for all time. When g k = (k € /C), no 
constraint is imposed on the real part, but the real-valuedness condition 13.251 requires 
that the imaginary component vanish. 
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Equation 13 . 371 can be solved by the method of integrating factors to obtain 

(3.38) u k (t) = ^/2D^pi f e~ ak( '- s) dB k ( S ). 

J — oc 

Since u k (i) is the projection of an Ito integral with deterministic integrand, it is at 
each time t a Gaussian random variable with mean zero. 
The variance of u k (i) can be computed from 13.381 by 

(3.39) E (|u k | 2 ) = Tr (e (u k fl^ 3 



To proceed further, we must distinguish between the cases k S JC and k ^ JC, where 
JC is defined in 13.291 

For k £ JC, we have pj^ = X and from constraint 13.361 that the Brownian motion 
Bk(s) is real valued. Therefore, 



(3.40) E [ dB k (s)dBl(s')) = 16(s - s')dsds', 
and it follows that 

(3.41) E (\urf) = ^ 
when k S JC. 

For k ^ JC, the Brownian motion is complex valued and pj^ is a projection onto 
the two-dimensional subspace orthogonal to gk- Therefore, 

(3.42) E (<2B k (s)dB£(s')) = 2p^S(s - s')dsds' 
and 

(3.43) Tr (p£) = 2, 
from which it follows that 



(3.44) S(|u k | 2 ) 



4Dk 
a k 



for k ^ JC. We remark that the formal calculations above can be justified rigorously 
by applying Ito's Isometry directly to equation 13.391 see reference |59l). 

To determine Z?k, we shall now compare these results with those that are obtained 
if we impose the condition that the immersed boundary method exhibit fluctuations 
governed by the Boltzmann distribution, as required by classical statistical mechanics. 
By Parseval's Lemma, the total kinetic energy can be expressed in terms of the Fourier 
modes of the fluid by 

(3.45) ^{uklH^lu^Ax 3 



k 
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The density of the Boltzmann distribution is then given by 
(3.46) *({u k }) = iexpf-^ 3 ^ |,ik|2 



Z \ 2k B T 

where Z is the partition function obtained by integrating Uk over the constrained 
subspace 



(3.47) n = {{u k } | g k • u k = 0, u N -k = u k }- 

Each degree of freedom of the fluid contributes a quadratic term to the energy of the 
system, giving a Boltzmann distribution which is Gaussian. Therefore, the equipar- 
tition theorem holds and each independent degree of freedom contributes on average 
^k B T to the kinetic energy. 

For a particular wavenumber kg£, the mean contribution to the energy is 

(3.48) 

where the expression 13.411 for E (|uk| 2 ) has been used. For such wavenumbers there 
are 3 independent degrees of freedom corresponding to the 3 real components of Uk- 
By the equipartition theorem this requires 

3.49 P - ^ = -k B T, 

2 ait 2 

which gives 

k B T 



pL 3 



(3.50) Ac = "k 

when k S JC. 

For k ^ AC, we must consider the pair (k, N — k) together, since fik = UN-k- The 
contribution to the mean energy of these two wavenumbers is 

(3.51, ^ 

where the expression 13.441 for E (|flk| 2 ) and E (|uN-k| 2 ) has been used. The number 
of independent degrees of freedom corresponding to the pair of wavenumbers (k, N— k) 
is 4, since the real vector space orthogonal to gk is two-dimensional and ilk is complex 
valued. 

The equipartition theorem in this case requires that 

(3-52) 2— = -k B T, 

2 oik 2 

which gives 

k B T 



(3.53) D k - a k 



2pL 3 



when k ^ JC. 
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To summarize, the following fluctuation-dissipation relation (|47t 1 701 ) is obtained 
when considering the constraints 13.251 and 13.221 imposed on the velocity field of the 
fluid: 



3 - 54 ^ k= ^» wr 

We remark that the Fourier mode of the fluid associated with k = [0,0, 0] T 
corresponds to translation of the fluid as a whole. From 13.231 the zero mode has 
«o = 0, which indicates that the fluid has no translational damping. As a consequence 
of 13. 541 the mode Uo is not thermally forced, which can also be understood physically 
by the conservation of total momentum by the internal thermal fluctuations. Thus 
for a fluid initially at rest with no net external force on the fluid as a whole, the 
translational mode remains zero Uo = under the thermal forcing. 

3.3.4. Numerical Method for the Fluid. To deal with the significant range 
in time scales for the modes of the fluid and immersed structures, we develop a time- 
stepping scheme that freezes the positions and forces exerted by the elementary par- 
ticles over a time step At, but otherwise integrates the dynamical equations exactly. 
With this approximation the set of equations 13.371 can be solved analytically using 
the methods of stochastic calculus (j59l ). This strategy has similarities to "exponential 
time differencing" or "exact linear part" numerical methods (3^; 36t l42j). 



In stochastic differential notation, the fluid equations with both thermal and 
particle forces can be expressed as 



(3.55) <iu k = -a k u k cft + p 1 p^f k dt + \/2L> k p k dB k (t), 

where to simplify the notation the subscript will be dropped for the Fourier modes of 
the particle force density so that f k always refers to fp r t, it- 
Approximating the particle force as constant over the time interval [t',t] gives 



(3.56) u k (i) = e-^-^u^t') + — (l - e-^t'A p^(t') 

+ VZD^J t/ e- ak( *- s) PkdB k ( S ), 

where p k is the projection operation defined in 13.281 and J •p k 'dB k (s) denotes inte- 
gration in the sense of Ito (|59l ) over the projected complex-valued Brownian motion 
Bk(£) defined in Section 13.3.31 

To obtain a numerical scheme for the fluid with finite time step At, each mode is 
updated at discrete times nAt using the analytic solution 13. 561 yielding the stochastic 
recurrence equation 



(3.57) u™ +1 = e-^'uH + — (1 - e-«* At ) p£f£ + Pk * 

po? k 



where uJJ = u k (nAt), f£ = f k (nAt), and H k = cr k J7 k . 

The notation fj k denotes a three dimensional complex-valued random variable, 
with each real and imaginary component being an independent Gaussian random 
variable with mean and variance 1. The random variable H k accounts for the 
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contributions of the stochastic integral in 13.561 over the time step. The variance <j\ 
can be determined by Ito's Isometry (|59l ) and is given by 

(3.58) al = ^ (1 - e - 2a - At ) . 

a k 

The constraint 13.251 that ensures the real-valuedness of the velocity field is re- 
spected by only applying the update T3. 571 to one member of each complex-conjugate 
pair, and then setting the new value for the partner mode as the complex conjugate of 
the computed mode. The condition 13.251 also requires that the modes with indices 
k £ K, have zero imaginary part; this is enforced explicitly in each time step. 

3.3.5. Numerical Method for the Immersed Structures. A time-discretization 
for the equation (|3.19p is developed for the advection of the elementary particles by 
integrating the fluid velocity field over a time step and then averaging the integrated 
velocity over a spatial neighborhood centered on the old particle position: 

(3.59) X" +1 'M - X">M = Vi„(x m - X"'M) f n+1 Um (s)dsAx 3 , 

where t n = nAt and X n 'M = X^(nAt). 

A precise integration of the fluid velocity u is taken which allows for time steps 
which underresolve the dynamics of some of the Fourier modes of the fluid. This 
capability is important due to the wide range of time scales that may be associated 
with the fluid modes and immersed structures in applications. If the time integral is 
approximated through numerical methods built from (stochastic) Taylor expansions 
about discrete times, such as Runge-Kutta methods and their stochastic variations 



(|44c l80l ). then it is important that the method sufficiently resolve the fluctuations 
of the processes to capture cancellations that occur over time. For instance, if the 
cancellation is not adequately captured, the numerical value of the integral of velocity 
will be larger in magnitude than the actual time integrated velocity. For immersed 
particles, this leads to an overly diffuse behavior where the particles overshoot their 
correct positions each time step. 

From 13.371 the time scale associated with the dynamics of the k" 1 mode of the 
fluid is l/«k- For the fastest modes of the fluid relevant for the immersed particle 
dynamics, the above considerations would place a severe restriction on the time step. 
While there may be clever numerical methods involving (stochastic) Taylor expansions 
which perform better than anticipated, a different approach will be taken here. 

To develop a method that remains accurate for a range of time steps, from those 
that fully resolve, partially resolve, or completely underresolve the fluid modes, we 
calculate the time integral in 13. 591 by substitution of the analytical expression 13 . 561 for 
the Fourier modes of the fluid velocity field. We recall that this approximation only 
assumes that the elementary particle positions and forces can be considered frozen 
over a time step. The resulting numerical scheme can therefore be expected to be 
accurate provided the time step At is chosen small compared to the time scales of 
the immersed structures, but with no restriction on the size of the time step relative 
to the time scales of the fluid modes. We will explain this property more precisely 
through numerical error analysis in Section |4] 

In updating the elementary particle positions in the numerical method, the time 



STOCHASTIC IMMERSED BOUNDARY METHOD 



17 



integral in 13.591 will be simulated as a random variable 

(3.60) = J " +1 u m (s)ds = ^ f JJ exp (i2vrk • m/N) 



(3.61) f £ = jf " +1 u k ( S )d S , 

with Ufe(s) given by (|3.56|) . 

By using standard techniques from stochastic calculus, the time integral can be 
evaluated by defining TJJ to give the Gaussian random variable 

(3.62) 




Using [3~56l at times nAi and (n + l)At, this can be expressed more simply as 



• / " +1 e-°*^-^p£dB k (r) + ^± ( Pk -B k (t n+1 ) - Pk "B k (i n ) 



(3.63) rg = -— (u»+! -u^) +p- 1 ^^Ai 



V2^ 



p k B k (i n +i) - p k B k (t„) 



The numerical scheme to update the elementary particle positions is then given 

by 

(3.64) x n+1 ^1 - X n 'W = <U x m - X n ^)r^Ax 3 , 

m 

where is generated each time step. To consistently update the particle positions 
with the velocity field of the fluid it is required that f 1 ^ be generated with the correct 
correlations to the fluid modes at the beginning and end of each time step, {u£} and 
{ u k +1 }- I n the next Subsection, a practical approach for doing so is presented. 

3.3.6. Method for Generating Modes of the Time Integrated Velocity 
Field. Since the modes TjJ of the time integrated velocity field and the modes u£ 
and u£ +1 of the velocity field evaluated at the beginning and end of a time step are 
not statistically independent, some care must be taken in generating the correspond- 
ing random variables that are used in the simulation. Since f^, u£, and u£ +1 are 
jointly Gaussian distributed random variables with mean zero, we need only ensure 
they have the correct covariances between their components. Since the real and imag- 
inary parts of each mode are independent, we shall for clarity consider only the real 
components with the understanding that the imaginary components are handled in a 
similar manner. 

In deriving a method to generate the time integrated field, it is useful to express 
the real part Re(r£) in terms of the following random variables 

(3.65) Re(f I) = p^Ao + p^A 1 + Pk "A 2 , 
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with 

1 _ p -a k At ( A f / i \ 2 \ 

(3.66) A = Re(u£) + — + — (e"^ At - l) ^Reflf) 

a k \ «k V ak / / 



(3.67) Ai = -Jl^Ji / e - Qk(t ™ +1 - s) Re(dB k (s)) 



c 

._Ls n 

a k k 



-.+1 



(3.68) A 2 = ^^Re ( B k (i n+1 ) - B k (t n ) 
a k 



+i 



Re(dB k (s)). 



The random variables were obtain by reorganizing the terms of 13.621 

This expression recasts the problem of determining the correlations of Re(rjJ) to 
the problem of determining the correlations of Ao, Ai and A2 with each other and 
the modes of the fluid velocity. A convenient feature of this approach is that Ao is 
already determined at the beginning of the time step, and is statistically independent 
of Ai and A2 by the independent increment property of Brownian motion. This 
reduces the problem to finding the covariance of Ai and A2. A useful identity for Ito 
integrals in this context is (|59l ) 

(3.69) e(J f(s)dB s J^ g{r)dB^j = J f(a)g{a)ds, 

where the notation E(-) denotes expectation with respect the underlying Brownian 
motion |59l). 

Using RT691 the covariance is given by 

(3.70) E(a[ j) A { 2 j) ) = -^(1- exp(-a k Ai)) , 

where the parenthesized superscripts denote the indices of the vector components. 
When j 7^ f the components and A^ are independent and have zero correla- 
tion. 

The variance of the components of Ai and A2 are given by 

(3.71) E(\A[ j) | 2 ) = ^ (1 - exp (-2a k Af)) 

(3.72) E (\A^f) = ^±At. 

From the numerical updating of the fluid variables described in Subsection l3.3.41 
Ai = — is already known each time step, so only A 2 need be generated. Ob- 

taining this random variable with the correct correlations can be accomplished by 
generating new standard Gaussian random variables ifi) (independent in j with mean 
and variance 1) and by taking the linear combination of the two random variables 
A^' and 77 ^ given by 



(3.73) 



A« =a 1 AP+a 2 7 7 « 
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with 



(3.74) 01 = v 



E[\A\ 



and 



(3.75) a 2 



E (|A^| 2 ) E (|A< j) | 2 ) - E (A^A^ 



In this manner, Re(f k ) can be generated from Ao, Ai, and A2 with proper 
accounting of correlations with the modes of the velocity field. The imaginary com- 
ponent Im(r k ) is generated in an analogous manner. 

4. Accuracy of the Method. In this section the accuracy of the numerical 
method is investigated. Three asymptotic scaling regimes of the time step are consid- 
ered. The first regime applies when the time step is taken sufficiently small to fully 
resolve the dynamics of the fluid. The second applies when the time step is taken 
large and completely underresolves the dynamics of the fluid. We finally consider the 
case in which the time step resolves some but not all of the fluid modes. 

Formal error estimates are given which show how the numerical errors scale with 
respect to the time step and various key parameters. While a rigorous analysis making 
use of standard stochastic Taylor expansion approaches (0) can be carried out for 
time steps which are small when compared to the time scales of the fluid and immersed 
structure dynamics, a completely rigorous analysis of the numerical error when the 
time step is large and underresolves a subset of the fluid modes is considerably more 
difficult. 

An important feature of the numerical method is the way in which the statisti- 
cal contributions of the fluid dynamics are taken into account, even when the fluid 
dynamics are underresolved. As discussed in Subsections 13 .3 .41 and [3~3 .51 the random 
increments of the elementary particle positions and fluid modes are simulated in such 
a way that the correct statistics and correlations are preserved over time steps which 
need only be small compared to the time scales of the immersed structures. While 
the time step relative to the time scale of the fastest modes of the fluid may be large, 
this procedure helps keep the local time discretization error small. By contrast, stan- 
dard finite difference schemes would generally have poor accuracy once the time step 
exceeded the time scales of the fastest fluid modes. 

To quantify the accuracy of the method, the strong error is considered, as defined 



in (44). Let X^l(t) denote the exact solution of equation 13.191 for the elementary 
particles and Uk(f) denote the exact solution to equation 13. 161 for the Fourier modes 
of the velocity field of the fluid. Let the numerically computed trajectories of the 
elementary particles be denoted by X^'] (t) an d the numerically computed fluid modes 
be denoted by Uk(t). Since we shall be interested in the error associated with a typical 
elementary particle, the superscript j will be dropped throughout the discussion. 

The strong error of the numerical method associated with the k th mode of the 
fluid is defined as 

(4.1) e M JAt) = E ( u k (At) - fi k (Ai) 
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and the strong error associated with an elementary particle is defined as 
(4.2) eprt(At) = E(|x(At)-X(At)|). 

The error associated to the velocity field in physical space is defined as 



(4.3) 



e M (At) =e(-^Y1 l u m(Ai) - iim(At)| Az 3 ^ 



For further discussion of the strong error see 0). 

The error expressions above and the estimates given below are intended to char- 
acterize the "typical" error for the numerical method; in reality they will of course 
depend on the particular configuration the elementary particles happen to be in at 
the beginning of a time step and the details of the forces acting between them. For 
the purposes of describing the errors incurred in the numerical method's handling of 
the force interactions, we shall therefore concern ourselves with describing how the 
errors scale with respect to the various numerical parameters. 

In the derivation of the estimates we quantify the error incurred by the numerical 
method's representation of the stochastic (thermal) components of the structural and 
fluid dynamics. The estimates presented follow from a systematic formal analysis of 
the errors resulting from the discretization of the stochastic and deterministic compo- 
nents of the dynamics, including their interaction during a time step. This calculation 
leads to a uniformly valid expression for time steps sufficiently small that the elemen- 
tary particles do not move appreciably (relative to their size) during a time step; no 
assumption is made in the derivation about the magnitude of the time step relative to 
the time scales of the fluid modes. As the resulting derivations are somewhat techni- 
cal, we shall in the present paper be content to state the error estimates, discuss their 
significance, and confirm their validity in a few special cases by numerical simulation. 
For a detailed derivation see jH). 

4.1. Error Estimates for Time Steps which Fully Resolve the Fluid 
Dynamics. When the time step is taken sufficiently small so that the dynamics 
of all modes of the fluid are resolved by the stochastic immersed boundary method 
(At <C min ^-), the following error estimates can be established: 

M F* / M 1 \ 

(4.4) g fld ,k(At) « —-Km (r + c n) ( c '^ + c'Vm) At 2 



(4.5) e fld (Ai) « (^ + C-) (C'ifee + C"v thm ) A? 



pa 3/2 L 3/2 a 



(4.6) e pTt (At) « (Qiv 2 thm + Cv bc v thm + Cvl c ) 



At 2 
a 



+ M ^(r + c "t;) + c""v thm) At\ 

pa \"f a J 

where M is the number of elementary particles, F* is the magnitude of the force 
acting on the elementary particles, and lp is the length scale associated with changes 
in the particle force of order F*. It will be assumed throughout that a < lp. The 
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factor 6* k is the magnitude of the Fourier coefficient for mode k of the function S a , 
averaged over all shifts (see Appendix[B|. We remark that 6* k « 1/L 3 for |k| <C L/a, 
while 5* k decays rapidly for |k| ^> L/a. 

In this notation the factors C which are superscripted with primes are approxi- 
mately independent of the physical parameters, and can be thought of as order unity 
constants. To avoid cumbersome notation and overly emphasizing the role of these 
factors the notation is reused in each equation, with the understanding that C denotes 
distinct factors for each estimate. The subscripted factors Q are also approximately 
independent of the physical parameters. They are distinguished since numerical val- 
ues will be estimated for these factors in order to make a comparison between the 
theoretical estimates and numerical simulations in the case that F* = 0. 

To simplify the expressions the following terms are defined t>thm = \/ ksT / pcfi 
and «f rc — F* / pa. The facto r Vt y,m can be interpreted via the equipartition theorem 
of statistical mechanics (47|; |7(J ) as the velocity scale of thermal fluctuations of an 
elementary particle of size a, since the associated mass will be proportional to pa . 
The term Vf rc can be interpreted as the velocity scale associated with the motion of 
a particle of size a in a viscous fluid when a force of magnitude F* is applied to the 
particle, since the friction coefficient of a particle is generally proportional to pa (flBh. 

The error estimates indicate that the stochastic immersed boundary method has 
strong first order accuracy as the time step is taken small. An error proportional 
to At 3 is included in ep r t(Ai) because its coefficient in certain circumstances can 
make it comparable to the At 2 . We remark that the reported proportionality of 
the errors with respect to M, the number of elementary particles, is based on a 
worst-case scenario where all M particles are clustered near each other. In general 
the error is expected to scale with a smaller factor reflecting the actual number of 
particles clustered in a region. Since this depends on details of the force interaction 
between particles, obtaining a more precise error estimate is technically involved and 
somewhat application dependent. While in practice the actual numerical error will 
likely be somewhat better than these factors indicate, we leave further refinements to 
future work in the context of specific applications. 

An important observation is that in the absence of forces on the immersed struc- 
tures (F* — 0), the fluid modes are simulated exactly (for the reasons discussed 
in Subsection 13.3.4)1 . Only the elementary particle dynamics incur a temporal dis- 
cretization error in this case (see Subsection I3.3.5[) , with the strong error incurred 
being of first order. A more conventional time stepping scheme based on finite dif- 
ferences would typically incur an error for the velocity mode iik which includes a 
contribution which scales as CW _3 / 2 ?; t j lm (akAi)™ +1 / 2 for some integer n. Such an 
error fails to remain small compared to the actual velocity change over a time step 
as soon as At > 1/au- For the numerical method developed in Section [3l the exact 
representation of the stochastic fluid dynamics, apart from the response to the forces 
exerted by the immersed structures, maintains better accuracy even as the time step 
underresolves the fluid dynamics. 

As demonstrated in Figure 14.11 the theoretical error estimate for the elemen- 
tary particles agrees well with numerical simulations in the absence of particle forces 
(F* = 0). The numerical results were obtained from simulations of the fluid-particle 
system with physical parameters in Table l4~2l In the comparison, the factors Q were 
computed from theoretical expressions arising in the derivation of the estimates, and 
their numerical values are given in Appendix [Dj It should be emphasized that the 
error estimates are stated as formal approximations, not as upper bounds. In the case 
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that there are particle forces, the discretization error depends on a number of details 
of the force structure, and therefore numerical comparison with simulations is left to 
future work in the context of specific applications. 



x 10 




0.2 0.4 0.6 0.8 

At (ns) 



Fig. 4.1. A comparison of the analytic error estimate for the elementary particle positions 
e P rf(At) given in equation \4-b\ with a numerical estimation of the error for the small to intermediate 
time step regime At <C parjji with F* = and parameter values given in Table \4~J%\ The dashed line 
denotes e prt (At). The data points denote the error as estimated from numerical simulations using 
the method. The error bars indicate one standard deviation of the sampled values. To compute 
the numerical estimates of the error, simulations were performed with each given time step At 
and compared with an ensemble of reference trajectories obtained from a high resolution simulation 
sufficient to resolve the dynamics of all N 3 modes of the fluid. The high resolution simulation had 
a time step of At = 10~ 3 ns < minl/aj,.) where mini/a^ = 3.94 X 10 -2 ns for the parameter values 
given in Table \4~H\ From the results, we see that the estimate given in equation ^. 6\ quantifies the 
error well for At <C pa 2 /fj, = 0.976ns. 



4.2. Error Estimates for Time Steps which Underresolve All Fluid 
Modes. We now present estimates for the error of the numerical method when the 
time step is taken large enough to underresolve all modes of the fluid, but always 
small enough to resolve the elementary particle dynamics: 

(4.7) max — < At < r mov (a). 

The notation r mov (a) denotes the time required for an elementary particle to move 
a displacement equal to its size a either by advection or diffusion. In this regime the 
following error estimates can be established: 

(4.8) 



STOCHASTIC IMMERSED BOUNDARY METHOD 



23 




(4.10) 



eprt(At) « Q 2 -At 



where D denotes the diffusion coefficient of an immersed particle (see Section 15. ip 
and the factors C and Q denote order unity nondimensional constants as discussed in 
Subsection 23] The other terms are the same as in Subsection 14. 11 

The smaller powers of At appearing in the error estimates may suggest that the 
accuracy is deteriorating more rapidly with respect to the time step in the underre- 
solved regime under discussion, but in fact the opposite is true. The error estimates 
reported above are in fact, for the range of time steps defining the underresolved 
regime, considerably smaller than the extrapolation of the error estimates in Subsec- 
tion [4J] which are valid only for the fully resolved regime. Indeed, the ratio of terms 
appearing in the above estimates to corresponding terms in the equations in Subsec- 
tion O involve ratios such as pL 2 / (p,\k\ 2 At) , D 1 / 2 /(w thm At 1 / 2 ) J and L x / 2 a?/ 2 p/ fj,At, 
all of which are much smaller than one in the asymptotic regime 14.71 

A more important point is that the numerical errors remain small relative to 
the changes in the system variables throughout this range of time steps, so that the 
numerical method maintains accuracy for all At < T mOY (a). This can be seen by 
observing that the changes in the system variables over a time step falling in the 
regime 14771 can be estimated as 

A/f F* A* 

(4.11) |M k |«C^^ + C'^§, 

(4.12) |H«C %C + C ; %m , 

(4.13) |(SX| « eVDAi + C'v frc At, 

where the factors C denote order unity nondimensional constants as discussed in 
Subsection 14. II The notation |<5[-]| indicates the absolute value of an increment of a 
variable over the time step. 

Since the velocity field of the fluid is completely underresolved, it changes by 
an amount comparable to its equilibrium value independently of the size of the time 
step. The ratios of the error estimates to the corresponding true changes in the 
system variables in 14. 131 can be bounded by sums and products of the nondimensional 
groups y/DAt/a, Vf TC At/a, M-^, and Mj^. The former two nondimensional groups 
involving the time step are both small by definition of the constraint At -C T mov (a) 
determining the asymptotic regime 14.71 under consideration. The nondimensional 
parameters M-^ and will be order unity or smaller when the system involves 

a small number of elementary particles. When the system contains a large number 
of elementary particles, these nondimensional groups can become large and the error 
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estimates become worse. While it is certainly to be expected that the presence of more 
complex structures involving more elementary particles will generally incur more error 
in the numerical simulation, we stress that the scaling of our errors with large M are 
surely too pessimistic. We therefore do not lay undue emphasis on the behavior of 
the errors for large M, which in any event will depend heavily on the details of the 
force structure. 

We emphasize that unlike traditional numerical analysis the presence of terms 
proportional to At in the error estimate 14.101 does not imply that the method is 
inconsistent. It must be remembered that these error estimates are appropriate not in 
the At I limit, but rather in the asymptotic regime [4~7l A more careful consideration 
of the sizes of the errors relative to the true changes in the system variables over a 
time step shows that our numerical method does in fact remain accurate for all time 
steps At <C T mov (a), even if the fluid modes are underresolved. Vital to this result 
was the use of the stochastic integral formula [3. 561 for the action of the thermal forces 
on the velocity held of the fluid, and the systematic consideration in Subsection l3.3.5l 
of how to correlate the stochastic component of the velocity field of the fluid with 
the random motion of the immersed structures. Without these developments, the 
resulting numerical method could not be expected to have good accuracy for time 
steps in the regime 14.71 

In Figure 14.21 the theoretical error estimate 14.101 for the elementary particle po- 
sitions over a long time step is compared with the results of a numerical simulation 
in the case that there are no particle forces (F* = 0). The numerical results were ob- 
tained from simulations of the fluid-particle system with physical parameters in Table 
14.21 In the comparison, the factors Q were computed from the theoretical analysis 
with values given in Appendix [D] The numerical simulations show good quantitative 
agreement with the formal error estimate 14.101 We remark that the estimate is to be 
understood as an approximation and not a rigorous upper bound. This agreement is 
evidence of the validity of the formal analysis. As discussed in Section [4~T| the errors 
arising in the presence of forces are not as explicitly quantifiable. We leave further 
discussion and verification of the estimates to future work in the context of specific 
applications. 

4.3. Error Estimates for Time Steps which Underresolve Only Some 
Fluid Modes. A key feature of the stochastic numerical scheme proposed in this 
work is that time steps can be chosen which only partially resolve the fluid dynamics. 
That is, the method need neither resolve all of the velocity modes nor comple tely 
neg lect the inertia of the velocity field (as in Brownian/Stokesian dynamics (fljl : I22I l75t 



l77l)). Rather, the time step can be chosen as needed to resolve the appropriate degrees 
of freedom of the fluid-particle system, having the fluid and thermal fluctuations 
interact appropriately with the structures. The case is now discussed in which the 
time step At falls within the intermediate regime 



where the dynamics of the fluid is only partially resolved. It turns out that the error 
estimates for the fully resolved regime (Subsection [41]) and the underresolved regime 
(Subsection 14. 2p each serve separately as formal upper bounds for all time steps, 
including the intermediate regime 14.141 Intuitively, then, one expects the numerical 
method to behave accurately over this intermediate range of time scales as well. To 
provide more quantitative support for this statement, error estimates are developed 




(4.14) 



mm 



— < At < max — 
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Fig. 4.2. A comparison of the analytic error estimate for the elementary particle positions 
e prt (At) given in equation \4-10\ with a numerical estimation of the error for the large time step 
regime maxl/ajj -C At -C r mov with F* = and parameter values given in Table \4~H\ The dashed 
line denotes e prt (At). The data points denote the error as estimated from numerical simulations 
using the method. The error bars indicate one standard deviation of the sampled values. To compute 
the numerical estimates of the error, simulations were performed with each given time step At and 
compared with an ensemble of reference trajectories obtained from a high resolution simulation with 
a time step of At = 1.0ns < maxl/a^ where maxl/a^ = 25.4ns for the parameter values given in 
Table \4^ From the results, we see that the estimate given in equation ^. 10\ quantifies the error well 
for max l/a^ <C At <C r mov , where T mov s;a 2 /D = 1.95 X 10 5 ns. 



for the asymptotic regime 

pa 2 „ 1 
(4.15) — <At<max — . 

it a k 

Then, for fluid modes that are well resolved, we have 
(4.16) 

%d,k(At) « ^p£,k (j£ + C~) (Cv bc At 2 + C'^DAt^) if a k At « 1, 

while the underresolved modes (with a^At ^> 1) have the same error estimate 14.81 as 
in the fully underresolved case. 

The errors incurred in the physical space variables describing the velocity and 
elementary particle positions can, in the asymptotic regime l4.151 be estimated as 



(4.17) 



(4.18) 
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Sp rt (At) « Q 2 —At 



a 
M 

ip ' a 



where v = p/p. 

These errors can be compared with the size of the actual changes in the system 
variables over a time step in the regime 14.151 which can be estimated by the same 
formulas as 14. 131 except that the resolved velocity modes have changes of approximate 
size 

(4.19) |*u k | « c"-^± At + C' W^ - 

The ratio of the errors to the corresponding magnitudes of the actual changes of the 
system variables over a time step is controlled by sums and products of the nondi- 
mensional quantities y/DAt/a, v hc At/a, (vAtf^/L 1 / 2 , Mf , and The former 
three remain small in the asymptotic regime [4. 151 under consideration, while the last 
two nondimensional groups (independent of time step) are related to our somewhat 
pessimistic bound on the force errors, as discussed in Subsection 14. II The numerical 
method is thereby shown to remain theoretically accurate within this intermediate 
asymptotic regime. In the absence of particle forces (F* — 0), the error estimates 
become identical to those for the unresolved fluid regime (Subsection l4.2p . 
One could also study the intermediate asymptotic regime 

1 . pa 2 

(4.20) min — < At < — , 

ak p 

which exists only when Ax <C a. For these time steps, all error estimates presented 
in Subsection 14.11 for the fully resolved regime remain valid, except that the estimate 
for the individual underresolved fluid modes is altered to 

A T F 1 * T 2 f A/f 1 \ 

(4.21) e MM (At) « — p-p- + C- J 8: M (CV C + C"v thm ) At. 

As with the other regimes, the errors in this regime are small relative to the magnitude 
of the changes of the actual system variables over a time step. 

By simple extension of the above arguments for time steps falling at the transitions 
between the asymptotic regimes, we see that the numerical method has been designed 
to remain theoretically accurate for all time steps At <C T mov (a), regardless of how 
well the fluid dynamics are resolved. 

5. Physical Behavior of the Method and Numerical Results. To ensure 
that the immersed boundary method with thermal fluctuations serves as a plausi- 
ble physical framework for modeling microscalc systems, we verify that the method 
exhibits several fundamental features which are correct according to the laws of sta- 
tistical physics |70l ) . In Subsection 15.11 an expression for the diffusion coefficient of 
immersed particles is derived, and it is shown that in three dimensions the mean 
squared displacement scales linearly in time and inversely in the particle size. It 
is further shown in Subsection 15.21 that the stochastic immersed boundary method 
captures the correct t -3 / 2 power law for the decay of the tail of the autocorrelation 
function of the particle velocity H; [H; 0; [H; [Hi; [H H3) ■ 
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In determining the thermal forcing in Subsection 13.3.31 we imposed the require- 
ment that the degrees of freedom of the fluid obey Boltzmann statistics in thermal 
equilibrium. In fact, the complete system including immersed structures should obey 
Boltzmann statistics. In Subsection 15. 3i we study the equilibrium statistics of im- 
mersed particles subject to a conservative force and show through numerical simu- 
lation that they do exhibit the correct Boltzmann statistics. To demonstrate some 
applications and as a further verification of the physical plausibility of the method, 
it is shown in Subsections 15.41 - 15.61 how the method can be used to model osmotic 



effects such as the pressure of confined particles, dimers, and polymers (|21t l6lh . In 
Subsection 15. 71 another application to a basic model of a molecular motor protein 
immersed in a fluid subjected to a hydrodynamic load force is presented (f66h . 

5.1. Diffusion of Immersed Particles. In this section the diffusion of particles 
in the stochastic immersed boundary method is discussed and an expression for the 
diffusion coefficient is derived. As part of the analysis it is shown that the correct 
diffusive scaling is obtained for three dimensional systems. To verify the validity of 
the approximations made in the analysis and to demonstrate the applicability of these 
results in practice, the results of the analysis are compared to the results of numerical 
simulations. 

In three dimensions, the diffusion coefficient for a single particle is defined as 

(5.1) D=lim (|X (1 )-X(0)|>) 

f->oo Qt 

In the notation the superscripts on the particle position are suppressed since only a 
single immersed particle is considered. 

An estimate for the diffusion coefficient of a single particle (with no interactions 
with other particles) in the stochastic immersed boundary method is derived in Sub- 
section [5TTT] from the autocorrelation function of the velocity field of the fluid. This 
estimate can be expressed as 

(5.2) j?=M^ E l^l 2 ^ 

3/3 ^ Ok 

where Yk is defined in appendix [Cl and (5 a ,k is defined in appendix [Bl This diffusivity 
as simulated by the stochastic immersed boundary method exhibits the physically 
correct scaling with respect to physical parameters (|45h . 

The diffusion coefficient is estimated from the numerical simulations using 

1 JL I 2 

(5.3) d^ — Y, |x m (*i)-X m (0) , 

m— 1 

where n is the number of sampled trajectories of fixed duration t\. The notation X m 
denotes the simulated particle position from the m th trajectory. 

In Figure 15.11 the theoretical estimate of the diffusion coefficient as given in 
equation l5.2l is compared to the numerical estimate given in equation [5?3] for particles 
with sizes a = 1, 2, 3, 4, 5 and for a long time step which underresolves the dynamics 
of the fluid. For the parameters of the fluid-particle system used in the numerical 
simulations, see table l4~2l For each particle size, the numerical estimates were made 
from n = 10 4 sampled trajectories with At — 10 3 ns and t\ = 10 4 ns. We remark that 
in the simulations, while the time step underresolves the fastest modes of the fluid, 
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Fig. 5.1. A comparison of the analytic estimate 1 5. 2\ for particles of the IB method with nu- 
merically estimated diffusion coefficients for the parameters given in Table \4-2\ The dashed line 
indicates the analytic estimate \5. 2\ The data points denote estimates of the diffusion coefficient ob- 
tained from numerical simulations using the method. The error bars denote one standard deviation 
of the sampled values. In the numerical simulations long time steps satisfying maxl/a^ <C At were 
taken which under For the parameters values given in Table \4^\ the time step At = 10 3 ns was used 
in the numerical simulations, where maxl/a^ = 25.4ns. The particle size ko corresponds to the 
parameter value a = koAx of the particle representation function S a defined in appendix HI where 
kg controls the number of mesh-widths spanned by the 8 a function. 



there is still good agreement between the diffusion coefficient of the simulated particle 
and the theoretical estimate 15.21 A tabulation of the theoretical diffusion coefficient 
as given in equation 15.21 for the immersed boundary method with various ratios of 
a I Ax and a/L can be found in (j45h . 

5.1.1. Derivation of the Diffusion Coefficient. To derive an analytic es- 
timate for the diffusion coefficient the semidiscrctizcd equations 13.161 and 13.181 are 
considered. Alternatively, the diffusion coefficient can also be derived directly from 
the stochastic immersed boundary equations by a stochastic mode reduction proce- 
dure (lih. In equations 13 . 161 and 13.181 the immersed particle dynamics are given by 

(5.4) ^) =U (X(t),i), 
where U is given in Fourier space by 

(5.5) U(x, t)=J2 £ 3 <l,k(x)u k (t) exp (z2/Tx • k/L) , 

k 

with the coefficient 5 a ,k(x) defined in appendix [Bl 

The autocorrelation function of the velocity of an immersed particle is 



(5.6) 



R(t, t + r) = (U (X(t), t) ■ U (X (t + r) , t + r)) 
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To estimate this, we use the fact that the statistics of U(x, f) are approximately 
shift invariant in both t and x, with the approximation improving as the spatial grid 
is refined. We further assume that the time scale associated with the fluid velocity 
is small relative to the time scale of the immersed particle motion, so that a particle 
moves a negligible distance relative to its size a during the correlation time of the 
fluid velocity. These approximations give 

(5.7) R(t, t + t) » (U(0, 0) • U(0, t)) := R{t). 

By applying (|C.5p . the autocorrelation function of the velocity of an immersed 
particle can then be expressed as 

(5.8) R{T) « L ^aJaM' (Uk(*) ' "k' (t + t)) 



k,k' 



V L e \SaM\ 2 3—e 



keA' 



El e -<*K\T\ 



V L 6 \5 aM \H^e-^\ 
k B TL 3 



]T|<kk| 2 T k e-^M, 



P k 

where Y k is defined in appendix [C] 

An important point for the numerical method developed in Section [3] is that this 
structure of the correlation function is preserved even for finite time steps, provided 
only that the time step is small enough that the immersed elementary particles do 
not move significantly during a time step (At <C T"diff(a)), where r diff (a) is the time 
scale of a particle to diffuse over a distance equal to its size a. Were we to have used 
instead a numerical method based on a stochastic Taylor expansion we would 
have to restrict the time step At to be small enough so that i?(r) is well approximated 
by a Taylor expansion for |r| < At, which would add the additional restriction that 
At <C 1/ak- Our more accurate representation for the fluid dynamics over a time 
step allows us to obviate this other condition, as demonstrated in Figure [5Tl 

A useful identity relating the autocorrelation function to the mean squared dis- 
placement of an immersed particle is 



(5.9) <|X(t)-X(0)| 2 ) = ^ 



ds Jn dr 

t 



= 2 1 R{r) ■ (t - r)dr. 



/o 



This allows for the diffusion coefficient to be estimated by the Kubo formula (|47f ): 



(5-10) D = IJ i?(r)dr ' 



By substituting the estimate 15.81 into 15.101 and evaluating the integral, the expression 
15.21 for the diffusion coefficient is obtained. 
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5.2. Algebraic Decay of Velocity Autocorrelation Function. For immersed 
particles diffusing in a viscous fluid, the particle motion is strongly coupled to the mo- 
tion of the fluid. As a particle moves along a particular direction, fluid is dragged 
along with it. When the particle changes direction, it is resisted by a viscous force 
arising from its motion relative to the nearby fluid with momentum related to the 
recent past of the particle's motion. This induces a somewhat stronger memory in 
the particle velocity than a standard model based on a constant Stokes drag would 
predict. In particular, a careful analysis of physical Brownian motion, including a 
more detailed model for the force between an immersed particle and the surrounding 
fluid, yields for D < p/p (0 13 [U H ^ S M) 



(5.11) R(r) « ^gg^r" 3 / 2 for r » pa 2 / p. 

The condition D -C pj p can readily be checked to hold for typical microbiological 
systems. 

For the stochastic immersed boundary method, it is shown in Subsection 15.2.11 
that this general behavior is recovered with 



(5.12) R(t) 



C m B ijTi t~ 3/2 for pa 2 /p < r < pL 2 /p, 



M 3/2 



where C IB = 47r 3/2 ■ The constant prefactors differ slightly due to the different ways 
particles are represented in the physical model and immersed boundary method. 

The restriction that r -C pL 2 /p for the t~ 3 / 2 scaling in the immersed boundary 
method is a finite size effect which should have an analogue for physical Brownian 
motion. For very long times where r 3> pL 2 / p the correlation function i?(r) decays 
exponentially, with rate governed by that of the lowest wavenumber modes in the 
Fourier series (|5.8p . However, by these times the autocorrelation function would 
already be very small so this very long time regime is of little practical interest. 

These results show that the decay of the particle velocity autocorrelation function 
in the stochastic immersed boundary method has the correct scaling with respect to 
time and physical parameters. 

5.2.1. Derivation of Algebraic Decay of Velocity Autocorrelation Func- 
tion. In this discussion, the reference to wavenumbers k implicitly indicates the 
value within the equivalence class of aliased wavenumbers such that each compo- 
nent |k( J )| < N/2. For this purpose one can choose any scheme to select a unique 
value when k lies on the boundary of this set. 

First observe from the scaling properties of Fourier transforms and the definition 
of S a from ([A3]) that 

(5.13) S aM sa 4,o = -^3 for |k| < L/a, 

and <5 aj k decays rapidly with respect to |k|o/X. Along with the fact that the high 
wavenumber components of R(t) decay at a faster rate ak than the low wavenumber 
components, it then follows from the Fourier series representation [B~8l for the particle 
velocity autocorrelation function that the sum will be dominated by the terms with 
|k| <L/a 

Over the intermediate asymptotic time interval indicated in I5.12[ the time t is 
small compared to the decay time l/a^ ~ pL 2 /p of the low wavenumber modes 
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|k| ~ 1, but large compared to the decay time l/a k ~ P a2 /p of the (relatively high) 
wavenumber modes |k| ~ L/a corresponding to the length scale of the particle. Com- 
bining these observations, there exists a time-dependent wavenumber scale k c (t) which 
satisfies 1 < k c (t) < N/2 such that e~ akt ps 1 for |k| < k c (t) and e -0!k W ps for 
wavenumbers such that |k| 3> k c (t). Consequently, over the intermediate asymptotic 
time interval, the Fourier series 15.81 is dominated by contributions from wavenumbers 
1 < ]kj < k c if) <C N/2. These observations allow us to make the following simplifying 
approximations over the time interval pa 2 / p <f< pL 2 /p: 

• The prefactors multiplying the exponential in each Fourier series term may 
be approximated by their low wavenumber limits: 

(5.14) |<5 a , k | 2 T k » 2/L 6 for |k| « N/2. 

• The decay rate in the exponential may be approximated for |k| <C N/2 by its 
low wavenumber asymptotics 

(5.15) a k re A|k| 2 ; A = 4tt 2 pp' 1 L~ 2 . 

• The Fourier sum may be extended to the full integer lattice, because with 
the replacement 15.151 the additional terms for large wavenumbers will be 
exponentially small and make a negligible contribution. 

• This Fourier sum over the integer lattice can be approximated by an integral 
over continuous k, because the dominant contribution comes from a large 
number of lattice sites 1 < |k| < k c (t), with k c (t) ^> 1. 

Applying these simplifications and then changing to spherical coordinates with 
radial variable k = |k|, we obtain 

(5.16) fl(r)« [ exp(-^|k| 2 r) dk 



3 P L 3 
nksT 



/>oo 

/ k 2 exp(-Ak 2 T) dk 
Jo 




pL 3 

_ 8irk B T 1 
~ pL 3 2 

where the second equality follows readily by using standard facts about Gaussians. 
In particular, the integral can be treated as the expectation of the second moment 
by introducing the standard normalization factor. Using 15.151 and simplifying the 
expression yields 15.121 

5.3. Equilibrium Statistics of Immersed Particles. For the particle-fluid 
system with the fixed temperature T, volume V, and number of elementary particles 
M, with the particles subject to a conservative force field, we have from statistical me- 
chanics that the equilibrium probability density ^ of the elementary particle positions 
should have Boltzmann statistics: 

where £ is the energy of a configuration of elementary particles. The factor Z is the 
normalization factor so that the density integrates to one. The Boltzmann distribution 
arises from the thermodynamic condition that the equilibrium probability distribution 
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Fig. 5.2. Boltzmann Distribution vs Numerical Equilibrium Distribution. The spherical po- 
tential energy has the parameters Ri = 125, R2 = 250, c = 6kgT / (R2 — Ri). The simulation 
was run with the parameters N = 16, L = 1000 nm, Ax = L/N, fi = 6.02 X 10 5 amu/(nm ■ ns), 
p = 602 amu/nm 3 , T = 300 K, At = 1000 ns. 200,000 time steps were simulated. 



of the microscopic states maximize entropy while maintaining a fixed average energy 
for the system (|70l ). 

In deriving the thermal forcing of the system in Subsection 13 . 3 .31 only the energy 
associated with the fluid modes was considered. When immersed structures are sub- 
ject to a force, it is not entirely clear that the correct equilibrium distribution for the 
system as a whole will be attained. 

When formulated in continuous time, the fluid-particle coupling in the immersed 
boundary method conserves energy exactly (64). While this may suggest that the 
correct equilibrium statistics should be obtained (up to the appropriate definition of 
an effective temperature) , the discretization of time in the numerical method could in 



rjrinciple disrupt it, particularly since we are not employing a symplectic method (|57 



74J). 

We now show that the numerical method from Section \3 . 1 1 appears to yield results 
consistent with the Boltzmann distribution, at least for the statistics of the position of 
a single immersed particle. For a more rigorous approach in which the Fokker-Plank 
equations associated with the stochastic immersed boundary method arc analyzed, 
see the related work (Q). 

To facilitate calculation of the equilibrium statistics, the particles are subject to 
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a radially symmetric external force depending on r = |X| with the potential energy 



It is assumed that c > so that this can be thought of physically as the potential 
associated with confining particles to a spherical chamber of radius R 2 - The inner 
radius R\ is used to soften the particle-wall interactions to avoid issues of numerical 
stiffness and R2 is taken significantly smaller than the spatial period of the lattice L. 
For a single immersed particle, the Boltzmann distribution for its radial coordinate 

is 



where Z is the normalization factor. 

In Figure 15.21 the equilibrium statistics of immersed particles simulated with the 
numerical method are compared with the Boltzmann distribution 1 5 . 1 9"! The simula- 
tions were performed with Ri = 125nm, R2 = 250nm and c = 6ksT/ (R 2 — Ri) with 
the parameters of the fluid-particle system given in Table 14.21 

5.4. Osmotic Pressure of Confined Non-interacting Particles. Osmosis 
is a phenomenon that occurs in many microscale biological systems. When diffusing 
particles are confined to a chamber by a boundary which is permeable to fluid but 
less permeable to particles, a pressure difference develops between the inside and the 
outside of the chamber. This difference is referred to as the "osmotic pressure" . 

When the confining boundary is impermeable to particles and the system is in 
equilibrium, van't Hoff's law (|70l ) relates the osmotic pressure to the concentration of 
the confined particles as 

(5.20) Posmosis = c k B T, 

where Co is the number of particles per unit volume in the chamber. More precisely, 
when the number of confined particles is small enough that the instantaneous pressure 
fluctuates, then van't Hoff's law should describe the ensemble or time average of the 
pressure difference that arises from confinement. 

One should see a signature of van't Hoff's law in the fluid pressure when a collec- 
tion of M non-interacting particles in a conservative force field with potential V are 
simulated by the stochastic immersed boundary method, given that the method was 
shown in Subection 15.31 to produce correct Boltzmann equilibrium statistics. Indeed, 
taking the expectation of the velocity, force, and pressure with respect to Boltzmann's 
distribution (ensemble average) in the fluid equation 12. II gives 



This is obtained using that (u) = and (f"thm) = 0- The notation (■) denotes 
the ensemble average over the thermal fluctuations and (p) denotes the average of the 
fluid pressure field. 

The ensemble average of the conservative force field with potential V at location 

x is 



(5.18) 




r < R x 

Ri < r < R 2 

r > R 2 . 



(5.19) 




(5.21) 



= -V(p(x)) +<f prt (x)>. 



(5.22) 




MW(x) 
Z 



k B T 
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Fig. 5.3. Time Averaged Radial Pressure. The spherical potential energy has the parameters 
Ri = 125, i?2 = 250, c = 6kgT/ (R2 — Ri). The simulation was run with the parameters N = 16, 
L = 1000 nm, Ax = L/N, fi = 6.02 X 10 5 amu/(nm ■ ns), p = 602 amu/nm 3 , T = 300 K, 
At = 1000 ns. 200, 000 time steps were simulated. 

The average pressure (p(x)) can be determined up to an additive constant from 
15.211 and 15.221 by taking the line integral in x. From the fundamental theorem of line 
integrals with the additive constant set to zero (to give the appropriate decay at large 
x) we have 

(5.23) (p(x)> = — e ^k B T 

By the assumption of Boltzmann statistics for the immersed particles the concentra- 
tion field is 

(5.24) c (x) = — e "bt. 

Substitution into 15.231 and integrating over the chamber gives van't Hoff's law. This 
derivation also indicates that the stochastic immersed boundary method has a fluid 
pressure field which should respect the local formulation of van't Hoff's law: 

(5.25) po(x) = c (x)k B T. 

From a statistical mechanical point of view of this derivation, this fluid pressure 
can be thought of as arising from the fluctuations of the system which persistently 
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subject a particle to the confining forces in the vicinity of the wall with a frequency 
determined by the Boltzmann statistics. The forces are then transferred to the fluid by 
the viscous particle-fluid interactions. Similar "mesoscopic" points of view of osmosis 
have been used in (|J; |2lj; 1611 ) . 

In Figure 15.31 a comparison is made of the pressure predicted by a local van't 
Hoff 's law taking into account the local solute concentration and the time average of 
the pressure field of the fluid obtained in a numerical simulation using the stochastic 
immersed boundary method for an immersed particle confined by the spherically 
symmetric potential given by equation 15. 181 

Another quantification of osmotic pressure is the average force per unit area 
exerted on the walls of the chamber which confine the solute. For a spherical chamber 
f2 of radius R in which the solute exerts a (generally repulsive) areal force density 
F pw (z) on a portion of the wall boundary <9f2 at relative location z, this osmotic 
pressure is given by: 

( 5 -26) Pwau = J F pw (y - x) • ^c(x) dy dx 

where c(x) denotes the average concentration of the solute particles. Note that this 
formula applies also when the solute particles interact with each other. We will 
only be considering isotropic wall-solute interactions (and uniform distribution of wall 
molecules) so that we can write the force density in terms of a (typically nonnegative) 
scalar function g pw (r): F pw (z) = q pw (|z|)z/|z| and the concentration density as c(x) = 
c(|x|). The expression in equation 15.261 for the osmotic pressure on the chamber wall 
can then be simplified by integrating over the angular degrees of freedom: 

(5.27) 2>wall = ^/ h(r)c(r)r 2 dr 
with 

fc(|x|)= f F pw (|y-x|).^dy 
Jan |y| 

7T r R+r 

(5.28) = - / q pw {p) (p 2 + R 2 - r 2 ) dp. 

r JR-r 

The function h{r) can be interpreted as the integrated normal force applied to the 
wall by a solute particle located a distance r = |x| from the origin. The change of 
variable used to obtain the last integral was p = \y — x|. 

When the potential confining the solute is "hard-walled," in the sense that the 
solute- wall interactions occur only in a very small boundary layer of the wall, the two 
formulas 15.271 and 15.231 give the same values for the osmotic pressure, up to a small 
difference which vanishes as the width of the boundary layer is taken to zero. For 
potentials which are " soft- walled" in the sense that solute molecules interact with the 
wall on a length scale comparable to the magnitude of the fluctuations of the size of the 
solute molecules, as in Subsection l5.5i the average pressure of the fluid and the average 
pressure exerted on the wall may in fact differ. In related work, we are investigating 
the various pressures associated with osmotic phenomena and exploring the influence 
of finite wall and molecule sizes (j!). For a discussion of how the osmotic pressure 
can be used to drive fluid flow in a mesoscopic pump, see the related work (|J). We 
now present a few examples to demonstrate how more complex structures immersed 
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in the fluid can be simulated with the stochastic immersed boundary method and to 
show how the osmotic pressure associated with wall forces can be derived from the 
thermal fluctuations of these structures as simulated by the method. 

Before proceeding in the following subsections to consider osmotic pressure effects 
of more complex structures, we remark that the wall pressure can be computed from 
the average of the radial confinement force f con f (x) = — / con f (|x|)x/|x| acting on the 
particles and equation 15.261 in the case of a spherical chamber. From Newton's third 
law (principle of equal and opposite forces): 

(5.29) fconf (x)= f -F pw (y-x)dy, 

Jon 

re-expressed using the isotropy of the forces involved, we have: 
(5-30) - /coBf(|x|)^7 = / -«(|y-x|)i^W 

i x i Jdo. \y ~ x i 

In spherical coordinates, this can be written: 

(5.31) / conf (r) = ^ [ R+7 q(p)(r 2 - R 2 + p 2 )dp. 

r JR-r 

If the model for the solute-wall interaction force can be assumed to vanish at 
separation distances comparable to the chamber radius, then this integral relation 
between q and / can be inverted to obtain: 



(5.32) 



\ L(R - P )f corit {R ~P) + (R + P)/co„f (R ~P) + J 9 fooni(s)ds 



\2irRp- 

The pressure on the wall can then be computed from the effective bulk confinement 
force /conf(^) using \5?27\ I5.28[ and this inversion formula. 

5.5. Application: Simulation of Interacting Immersed Particles and Os- 
motic Pressure. We now discuss application of the stochastic immersed boundary 
method in determining the osmotic pressure when the confined particles can inter- 
act. In particular, we consider the case in which particles interact in distinct pairs 
(dimers) through a spring with non-zero rest length and are confined to an approx- 
imately 400nm spherical chamber. Note that the solute particles are confined in 
a microscopic chamber, in the sense that the chamber diameter is comparable or 
smaller than the length-scale associated with the solute particle interactions between 
the monomers. This is in contrast to a macroscopic chamber in which solute particles 
interact on a length-scale very much smaller than the chamber diameter and where 
the van't Hoff law is well established with the osmotic pressure depending only on the 
number of solute particles and not on their physical characteristics. For example, in 
microscopic chambers the amplitude of the fluctuations of a solute particle's diameter 
may be comparable to the chamber diameter and play a non-negligible role in the 
osmotic pressure associated with confinement. 

To investigate these effects, we consider how the osmotic pressure changes as the 
binding strength for a collection of dimers is varied. From the classical van't Hoff 's 
law, it would be expected that the osmotic pressure for tightly bound dimers are half 
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that of the zero-binding case (free monomers) , because a tightly bound dimer behaves 
effectively as a single particle. This concept is exploited to suggest the design of an 
osmotically driven pumping apparatus in (|J). From simulations using the stochas- 
tic immersed boundary method, we can investigate how the osmotic pressure varies 
between the unbound and tightly bound regimes. 




Fig. 5.4. Illustration of five distinct pairs of coupled particles confined in a spherical chamber 
with the soft-wall potential given by equation \ 5.34\ Simulations were performed for five pairs of 
particles coupled with interaction energy given by equation 1 5. 3 'A 




none 



250 300 350 400 

r(nm] light 



Fig. 5.5. The geometrically weighted particle concentration c(|x| = r)r 2 of the monomers in 
the spherical shell of radius r used in equation \ 5.27\ As the coupling is tightened, the particle 
concentration has more weight at smaller radii. The concentration is shown only over the boundary 
layer over which the confinement force is exerted. 

Each pair of particles is coupled by a potential energy corresponding to a standard 
spring model with finite rest length I: 

(5.33) $ 2 (X 1 ,X 2 ) = |(|X 1 -X 2 |-^) 2 , 

where Xi and X 2 denote the particle locations and K represents the spring stiffness. 
Each monomer at location x is subject to a confinement force given by the radially 
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symmetric potential: 



(5.34) 



<&i(x) 



0, 

Qa 

4 (^2 - Rl 



x 



Ri 



\*\<Ri, 

Ri < |x| < R 2 

|x| > R 2 . 



This potential can be thought of as arising from the interaction force of the confined 
solute monomers with particles distributed uniformly over the walls of a spherical 
chamber having radius R — R 2 . The formula given in equation 15.321 gives the re- 
lationship between the monomer-wall interaction force and the radial confinement 
force. For the simulations the parameters were chosen as I = lOOnm, R± — 375nm, 
R 2 = 400nm, C = 8k B T/(R 2 - Ri) 2 ■ The setup is depicted pictorially in Figure l5Tl 
and movies of the simulations can be found in the Supplemental Materials. 
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Fig. 5.6. The osmotic pressure for confined dimers as a function of coupling strength. As the 
coupling strength of the dimers increases, the osmotic pressure decreases in a non-linear manner. 
The osmotic pressure value for unbound monomers is approximately double that of the the tightly 
bound monomers, which is in accordance with what is expected under van't Hoff's law. Deviations 
from van't Hoff's law are apparent for intermediate coupling strengths for which the length scale of 
the dimers is comparable to those of the chamber and the wall thickness. 



Figure I5T5I displays the results of simulations with the stochastic immersed bound- 
ary method that show that as the stiffness is increased, the particle density decreases 
for each radius r within the region of the confining potential. The pressure conse- 
quently drops with increasing coupling stiffness, as shown in Figure [5TB1 We see also 
that, in accordance with van't Hoff's law, the pressure in the strong coupling limit, 
where the particle pairs behave effectively as single entities, is cut to roughtly half 
from the no coupling case. We observe deviations from the classical van't Hoff law 
when the length scale of the bound molecules is comparable to that of the wall or the 
chamber. 

An intuitive statistical explanation for the pressure drop is that a strongly cou- 
pled particle is less likely to venture far into the confining potential because roughly 
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half the time its partner will encounter the confining potential first, be repelled, and 
pull its accompanying particle back away from the confining potential sooner than it 
would have on its own. In more physical terms, the entropy of the particle pairs de- 
creases and consequently the entropic penalty associated with confinement is reduced 
as the coupling strength increases. For coupling values that make the dimer length 
scale comparable to the microscopic chamber size, the osmotic pressure assumes an 
intermediate value which is not well described by a van't Hoff's law. 

5.6. Application: Simulation of Polymer Chains and Polymer Knots. A 

fundamental feature of the stochastic immersed boundary method is that each struc- 
ture evolves according to a local average of a common fluid velocity field. The method 
therefore automatically captures the physical phenomenon that the velocities of im- 
mersed structures become strongly correlated when they are close together in space. 
Mathematically, the solution map of the immersed structures and surrounding fluid 
volume, which maps a configuration of the fluid and structures at a reference time 
to the solution configuration at a later time t, can be viewed as a homeomorphism. 
Consequently, in the continuous-time framework, the method preserves topological 
invariants of the immersed structures, such as the knottedness of a continuous closed 
curve, as they evolve. This is in contrast to other simulation methods, such as Stoke- 



sian Dynamics (|15c l75c 1771 ) , which would require explicit excluded volume constraints 



and/or repulsion forces between monomers to prevent topological changes. 




230 2*0 260 280 270 2S0 2S0 300 
r(nm) 



Fig. 5.7. The bar graph shows the average concentration of monomers within the chamber 
in the boundary layer in which they interact with the confinement forces (Subsection \5.4P - As 
the knottedness of the polymer increases, the monomers are more restricted and concentrate on 
average toward the chamber center spending less time interacting with the confinement forces. As a 
consequence the concentration of monomers in the boundary layer decreases and the overall pressure 
drops (Table\5^Tj) . The parameters in the simulations were taken the same as in Subsection \5.4\ 



To demonstrate this feature of the method in practice (with finite time step) and 
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to show how worm-like chain polymers can be simulated, the stochastic immersed 
boundary method was applied for a generic polymer chain, a polymer trefoil knot, 
and a polymer figure eight knot. From these simulations the osmotic pressure of 
confinement was estimated for each of the polymers. The results of the average 
concentration of the polymer monomers, which determine the average radial force 
density exerted on the confining wall, is given in Figure 15.71 

As the knottedness of the polymer increases, its constituent monomers spend less 
time at large radii, and as seen in Table 15.11 the osmotic pressure is significantly re- 
duced. An intuitive explanation is that as the knottedness of the polymer increases, 
this restricts the intrinsic configurations accessible to the thermally fluctuating poly- 
mer. In physical terms, the knottedness reduces the entropic penalty of confining of 
the polymer. Movies showing simulations of the thermally fluctuating polymer knots 
can be found in the Supplemental Materials. 

Table 5.1 
Osmotic Pressure of Polymer Knots 



Knot Type 


Osmotic Pressure (amu/nm • ns 2 ) 


Unknotted 
Trefoil Knot 
Figure Eight Knot 


0.16 

0.0439 

0.0392 



5.7. Application: Simulation of a Basic Model for a Molecular Motor 
Protein Transporting a Membrane-Bound Cargo Vesicle . We now discuss 
how more complex systems can be simulated with the immersed boundary method. 
On a subcellular level motor proteins interact with cytoskeletal structures, such as 
actin and microtubules, to generate force and to transport materials within the cell. 
For example, neurotransmitters are produced in the cell body of neurons and trans- 
ported by kinesin motor proteins along axons to the vicinity of the synaptic cleft 
where they are packaged for future release Q). We demonstrate how the stochastic 
immersed boundary method can be applied to simulate a basic model of a molecular 
motor protein immersed in a fluid moving along a filament which transports a cargo 
vesicle (Figure I5~5| . 




Fig. 5.8. Illustration of the basic model of a molecular motor protein immersed in a fluid and 
towing a cargo vesicle. As the fluid flow strengthens, the hydrodynamic drag on the cargo increases 
and a load force is exerted in opposition to the motor transport. For large opposing fluid flows, the 
cargo may significantly change shape in response to the flow. 
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The motor protein is modeled as a Brownian Ratchet (41; 66) and the cargo vesicle 
is modeled by a triangulated mesh which forms a membrane enclosing a spherical 
volume. The nodes of the mesh are linked together by springs of the form given in 
equation 15.331 with non-zero rest lengths determined by the distance between nodes 
in an initial spherical configuration. The cargo is linked at the vesicle surface to the 
motor by a spring of the form given by equation 15.331 with a non-zero rest length of 
approximately 100 nm. The spherical vesicle has a radius of 125nm and the ratcheting 
intervals (light and dark inset of Figure I53|) are of length 100 nm. 




0.1 015 
drag force fpNi 



Fig. 5.9. The mean motor velocity vs the hydrodynamic drag force. The data points joined by 
the solid curve shows the mean motor velocity X t /t obtained from simulations over approximately 
t = 3ms. The dashed curve shows the mean velocity for an idealized ratchet having a viscous drag 
comparable to the spherical vesi cle and subject to a constant load force the same strength as the 
hydrodynamic drag force Liil . lfid) . 



To examine how the mean velocity of transport by the motor protein behaves un- 
der different loading conditions, simulations were performed in which a hydrodynamic 
drag force is generated on the vesicle cargo by a bulk flow of the fluid (@;El). The mean 
velocity for different strengths of the countering fluid flow is plotted in Figure I5.91 
where we see that for significantly large opposing flow, the motor can almost be made 
to "stall" . The stochastic immersed boundary method allows for hydrodynamic effects 
associated with the shape and deformation of cargos to be investigated. These effects, 
not typically considered in other numerical simulation approaches for molecular mo- 
tors, may have important consequences for motor/polymerization ratchet transport 
when more details of the motor are taken into account or with more complex cargos 
such as membrane tubes, small cell organelles, or chromosomes (@; f60h . More 
sophisticated models can also be formulated within the stochastic immersed bound- 
ary method framework, such as the case in which multiple motor proteins transport 
a common cargo, interact by crosslink cytoskeletal filaments, or include additional 
mechanical degrees of freedom of the motor protein itself. Some movies of our motor 
protein model transport simulations can be found in the Supplemental Materials. 
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6. Conclusion and Discussion. In this work we have discussed how thermal 
fluctuations can be incorporated into the immersed boundary method in a manner 
consistent with the laws of statistical mechanics. A new stochastic numerical method 
was proposed that allows for a range of time steps to be taken in which the fastest 
degrees of freedom of the fluid-particle system are either undcrresolved, partially 
resolved, or fully resolved. In addition, the numerical method was designed to take into 
account in a systematic way the statistical contributions of the thermal fluctuations 
over long time steps with the correct correlations between the particles and fluid. 

To investigate the behavior of the immersed boundary framework and the stochas- 
tic numerical method with respect to well-known laws in statistical physics a number 
of theoretical results were obtained for the method and compared with numerical 
simulations. In particular, it was shown that immersed particles simulated with the 
numerical method exhibit the correct scaling in the physical parameters for the mean 
squared displacement in three dimensions. It was also shown that the stochastic nu- 
merical method captured inertial effects of the fluid with a velocity autocorrelation 
function for a particle that for long times decays with algebraic order r~ 3 / 2 . We 
further found that particles appear to have the correct Boltzmann equilibrium statis- 
tics. Moreover, the method was found to produce the van't Hoff law of osmosis for 
a particle confined to a spherical chamber recovering the correct osmotic pressure. 
In addition results were presented which showed how the osmotic pressure could be 
computed for interacting pairs of particles and worm-like chain polymers, including a 
trefoil and figure eight polymer knot. A more complex application of a basic model 
of a molecular motor protein immersed in a fluid and towing a vesicle-bound cargo 
subject to a hydrodynamic drag force was simulated and the force- velocity statistics 
computed. 

These basic physical checks indicate that the stochastic immersed boundary method 
has the capability to capture many important features of thermally fluctuating sys- 
tems involving immersed structures which interact with a fluid. The results presented 
suggest the promise of the stochastic immersed boundary method as an effective ap- 
proach in modeling and simulating the mechanics of biological systems at the cellular 
and intracellular level. 
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Appendix A. The Representation Function 5 a for Immersed Particles. 

In the immersed boundary method, it is required that a function <5 a be specified to 
represent the elementary particles. The representation of this function is often derived 
from the following function which is known to have desirable numerical properties 
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(A.l) 



i (5 + 2r - V-7- 12r-4r 2 ) 



i (3 + 2r + s/l - 4r - 4r 2 ) 



± (3 - 2r + s/l + 4r - 4r 2 ) 



, if r < -2 
, if -2 < r < -1 
, if -1 < r < 
, if < r < 1 



|(5 - 2r - V-7 + 12r - 4r 2 ) , if 1 < r < 2 



I o 



if 2 < r. 



For three dimensional systems the function S a representing elementary particles 
of size a is 



(A.2) 



1 



,(i) 



,(2) 



.(3) 



where the superscript indicates the index of the vector component. 

To maintain good numerical properties, the particles are restricted to sizes a = 
nAx, where n is a positive integer. For a derivation and a detailed discussion of the 
properties of these functions see |64T). 

Appendix B. The Fourier Coefficients of the Function S a Used to Rep- 
resent an Immersed Particle. Throughout the paper it will be useful to consider 
the Fourier coefficients of the function <5 a (x — X) used to represent an elementary 
particle situated at position X. While the function is defined for all x € A, it is 
often useful to consider the restriction of the function to the discrete lattice points 
{x m = mAir|m g 

We will use the following notation to denote the discrete Fourier transform of the 
delta function restricted to the lattice: 



(B.l) 



<5a,k(X) = ^X^ a ( Xn 



X) exp (-i27rk ■ m/AT) . 



The dependence of the Fourier coefficients on the particle position X (relative to the 
lattice) is explicitly noted. When the dependence on X is not explicitly noted, then 
we will be referring implicitly to the discrete Fourier transform of the delta function 
when centered on a lattice point: S a ,k '■= S a ,k(0). 

Appendix C. Autocorrelation Function for the Velocity Field of the 

Fluid. In this section the autocorrelation function is computed for the velocity field 
of the fluid in the absence of force f pr t = 0. This is done by representing the velocity 
field in Fourier space and computing the autocorrelation function of each mode k. 
From equation 13. 561 and standard stochastic calculus the steady-state autocorrelation 
function of the k th mode when s > r is 




E (u k (s) ■ u k (r) 



e -a k (r-w) 



J —OO J 
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where the notation pj^ denotes projection orthogonal to gk as defined in Subsection 

MM 

By applying Ito's Isometry to IC.21 and observing the symmetry under the inter- 
change s <-> r, the autocorrelation function is given by 

(C3) E (u k ( S ) • u k (r)) = 4 ^ e _ ak | s _ r , 



where 
(C.4) 




3, ke /C 
2, k££, 



and the index set /C is defined in 13.291 

The factor Tk arises from the incompressibility constraint l3.22[ the real-valuedness 
constraint 13.251 and the dimensionality of the space orthogonal to g k . See Subsection 
13.3.31 for a discussion of how the constraints affect -D k . 

The spatio-temporal correlation function of the velocity field u is then given by 



E (u m (s) ■ u n (r)) = J2J2 E ( flk ( s ) ' Uk '( r )J cxp ( i27r ( n ' k ' m ' k )/ N ) 

k k' 

(C.5) = ^ E (u k (s) • u k (r)) exp (i27r(n - m) • k/iV) 



k 



P L3 k 



To obtain the second equality, we used the statistical independence of the Fourier 
modes of the velocity field when the indices k and k' are distinct and do not correspond 
to conjugate modes (see 13.25]) . When the indices k and k' do refer to conjugate but 
distinct modes, then the average vanishes because a mean zero random variable Z 
with independent and identically distributed real and imaginary components satisfies 
(Z 2 ) = 0. The last equality follows by substitution from equation IC.31 



Appendix D. Constants: Accuracy and Error Estimates. The non- 
dimensional factors Q appearing in the error estimates in Section [4] are approximately 
independent of the physical parameters. For comparison of the theoretical error es- 
timates with numerical simulations, it is useful to compute the factors for specific 
physical parameters to obtain estimated values. Evaluation of the expressions for the 
Q constants for the system with parameters given in Table 14.21 gives the following 
values: 

(D.l) Qi = 0.563 

(D.2) Q 2 = 7.87. 
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Table 4.1 
Parameters of the Method 



Parameter 


Description 


k B 


Boltzmann's constant 


T 


Temperature 


L 


Period Length of Fluid Domain 


P 


Fluid Dynamic Viscosity 


P 


Fluid Density 


N 


Number of Grid Points in each Dimension 


At 


Time Step 


Ax 


Space Between Grid Points 


a 


Effective Elementary Particle Size (approximate radius) 



Table 4.2 
Values used in Numerical Simulations 



Parameter 


Description 


T 


300 K 


L 


1000 nm 




6.02 x 10 5 amu/(nm • ns) 


P 


602 amu/nm 3 


N 


32 
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Table 4.3 
Notation Conventions 



Parameter 



Description 



6a 
5a,k 

a k 
D k 

U m 
Ufe 

U 

Xbl 

fprt 

f thm 

ffc 



Representation function of an immersed elementary particle of size a 
The k" 1 Fourier coefficient of the particle representation function 
Damping of the k. th Fourier mode 

Strength of the thermal forcing of the k. th Fourier mode 

Fluid velocity at the m th grid point 

The k*' 1 Fourier mode of the fluid velocity field 

Smoothed fluid velocity field for immersed elementary particles 

Position vector of the m th Eulerian grid point 

Position vector of the j th immersed elementary particle 

Force density arising from the immersed structures 

Force density arising from the thermal forcing 

The k th Fourier mode of the structural force density field 



